2.1 Starting with output from deepTools ‘computeMatrix’

If you do do not currently have deepTools installed, instructions for installation are in the deepTools manual, with the easiest method likely being through Bioconda.

This direct output from the ‘computeMatrix’ function can be imported into R as a profileplyr object using the import_deepToolsMat() function. This function takes as its only argument the path to the matrix from deepTools. The information contained within this matrix is stored in a profileplyr object.

library(profileplyr)

example <- system.file("extdata", 
                       "example_deepTools_MAT", 
                       package = "profileplyr") 
proplyrObject <- import_deepToolsMat(example) 

## [1] "Matrix contained 'NA' values, these will be replaced with zeros"

2.2 Starting with output from soGGi

The soGGi function plotRegion() allows for the quantification of signal over genomic intervals within R, and the information from this quantification is stored in a ChIPprofile object. profileplyr can take the ChIPprofile object as an input to generate a profileplyr object using the as_profileplyr() function. While the plotRegion function from soGGi does not allow for inputting multiple signal files or multiple BED files, profileplyr contains a function to do this, BamBigwig_to_ChIPprofile(). This function takes a character vector of paths to bigWig files or BAM files in the ‘signalFiles’ argument. Note that the files for each function call must be all BAMs or all bigWigs and not a combination of the two. The ‘testRanges’ argument is must be a character vector of paths to BED files. The corresponding names for each BED file that will appear on visualizations involving the resulting profileplyr object can be set with the ‘testRanges_names’ argument.

library(soGGi)

signalFiles <- c(system.file("extdata", 
                             "hindbrain1_binned.bw", 
                             package = "profileplyr"), 
                 system.file("extdata", 
                             "liver1_binned.bw", 
                             package = "profileplyr"),
                 system.file("extdata", 
                             "kidney1_binned.bw", 
                             package = "profileplyr"))

testRanges <- system.file("extdata", 
                          "ranges.bed", 
                          package = "profileplyr")

chipProfile <- BamBigwig_to_chipProfile(signalFiles, 
                                        testRanges, 
                                        format = "bigwig")
chipProfile

## class: ChIPprofile 
## dim: 763 3001 
## metadata(3): names AlignedReadsInBam group_boundaries
## assays(3): '' '' ''
## rownames(763): giID197 giID198 ... giID762 giID763
## rowData names(4): name score giID sgGroup
## colnames(3001): Point_Centre-1500 Point_Centre-1499 ...
##   Point_Centre1499 Point_Centre1500
## colData names(0):

Once we have the ChIPprofile object, the as_profileplyr() function will generate the profileplyr object that can be used in the many downstream functions described below. Similar to when you import from deepTools, the grouping of the ranges after using the as_profileplyr() function is stored in the ‘rowGroupsInUse’ section of metadata(proplyrObject). Further, if the object is derived from a soGGi ChIPprofile object, the inherited group column will be called ‘sgGroup’.

proplyrObject <- as_profileplyr(chipProfile)
metadata(proplyrObject)$rowGroupsInUse

## [1] "sgGroup"

While the as_profileplyr() function provides a convenient way to go straight from BAM and bigWig files to profileplyr objects, and allows the user to use any argument that can be passed to the soGGi regionPlot() function, it is possible to directly use the regionPlot() function. This can be done in one line of code using pipes (shown below), however, note that the as_profileplyr() function adds the ability to use multiple BED files or GRanges with one command, which cannot be accomplished in the code below.

library(BiocParallel)

testRanges <- import.bed(system.file("extdata", 
                                 "ranges.bed", 
                                 package = "profileplyr"))

proplyrObject <- bplapply(signalFiles,
                          regionPlot,
                          testRanges = testRanges, 
                          format = "bigwig",
                          style="percentOfRegion") %>% 
  do.call(c,.) %>% 
  as_profileplyr(names = "giID")

## [1] "Matrix contained 'NA' values, these will be replaced with zeros"

3.1 Key components of the object

The matrices that represent the bigWig or BAM signal within each bin (bin size specified within deepTools, soGGi by default uses base pair resolution, so bin = 1) are contained in a list within the profileplyr object that can be accessed with the assays() function. The columns of each matrix are the bins, while the rows are the ranges from the BED files that were used as the input to deepTools or soGGi, and the dimensions of each matrix can be seen with dim().

library(SummarizedExperiment)

assays(proplyrObject)

## List of length 3
## names(3): hindbrain1_binned liver1_binned kidney1_binned

dim(proplyrObject)

## [1] 763 100

A key feature of the RangedSummarizedExperiment class is the ability to link the rows of the assay matrices to genomic ranges. Those ranges are stored in a GRanges object allowing for efficient integration with other Bioconductor packages. The GRanges object within the profileplyr object can be accessed with the rowRanges() function. Standard GRanges accessor functions such as start(), end(), or ranges() can be used on the entire profileplyr object. Furthermore, as these ranges are annotated by the various functions of profileplyr, information and classifications of these ranges will be stored in this GRanges object in the metadata columns, which can be accessed as a Dataframe with the mcols() function (this same information is also contained in the rowData section of the object, accessed with rowData()).

rowRanges(proplyrObject)[1:3]

## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames            ranges strand |    names     score  dpGroup
##          <Rle>         <IRanges>  <Rle> | <factor> <numeric> <factor>
##   [1]     chr8 33902396-33903372      * |        .         0    genes
##   [2]    chr17 27678695-27679778      * |     ._r1         0    genes
##   [3]    chr17 23546216-23547490      * |     ._r2         0    genes
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths

The information associated with each sample is stored as a Dataframe within with the metadata of the profileplyr object, and is accessed with sampleData().

sampleData(proplyrObject)

## DataFrame with 3 rows and 19 columns
##                    upstream downstream      body  bin.size unscaled.5.prime
##                   <numeric>  <numeric> <numeric> <numeric>        <numeric>
## hindbrain1_binned         0          0      1000        10                0
## liver1_binned             0          0      1000        10                0
## kidney1_binned            0          0      1000        10                0
##                   unscaled.3.prime     sample_labels   verbose bin.avg.type
##                          <numeric>          <factor> <logical>     <factor>
## hindbrain1_binned                0 hindbrain1_binned     FALSE         mean
## liver1_binned                    0     liver1_binned     FALSE         mean
## kidney1_binned                   0    kidney1_binned     FALSE         mean
##                   missing.data.as.zero     scale skip.zeros nan.after.end
##                              <logical> <numeric>  <logical>     <logical>
## hindbrain1_binned                FALSE         1      FALSE         FALSE
## liver1_binned                    FALSE         1      FALSE         FALSE
## kidney1_binned                   FALSE         1      FALSE         FALSE
##                   proc.number sort.regions sort.using ref.point
##                     <numeric>     <factor>   <factor>    <list>
## hindbrain1_binned           1         keep       mean      NULL
## liver1_binned               1         keep       mean      NULL
## kidney1_binned              1         keep       mean      NULL
##                   min.threshold max.threshold
##                          <list>        <list>
## hindbrain1_binned          NULL          NULL
## liver1_binned              NULL          NULL
## kidney1_binned             NULL          NULL

The ‘rowGroupsInUse’ accessor for the metadata indicates which column of the range metadata (mcols() or rowRanges()) that will be used for grouping in the final output if that object were used for visualization. For a profileplyr object created from a deepTools computeMatrix output or from a soGGi ChIPprofile object, the inherited groups correspond to the BED files over which the signal was quantified, and these groups are contained in the ‘dpGroup’ column.

metadata(proplyrObject)$rowGroupsInUse

## [1] "dpGroup"

The ‘mcolToOrderBy’ slot of the metadata indicates which column of the range metadata will be used for ordering the ranges as they are exported to either deepTools or to EnrichedHeatmap (with the generateEnrichedHeatmap() function within profileplyr). This can be set using the orderBy() function, which requires a profileplyr object and a character string matching a column name of the range metadata as arguments. If groupBy is never used on an object, the default will be to order by the mean signal of each range (within each group). In addition, until groupBy() has been used on a profileplyr object, the ‘mcolToOrderBy’ slot of the metadata will be NULL and will not be seen with metadata(proplyrObject).

NOTE: If you are exporting to deepTools and want to order by something other than mean range signal (or other deepTools ordering options), then make sure to set the –sortRegions argument of deepTools ‘plotHeatmap’ to “no” so your custom ordering will be used. generateEnrichedHeatmap() will always use the column indicated by metadata(proplyrObject)$mcolToOrderBy, or the mean signal if that value is NULL.

metadata(proplyrObject)$mcolToOrderBy

## NULL

proplyrObject <- orderBy(proplyrObject, "score")
metadata(proplyrObject)$mcolToOrderBy

## [1] "score"

3.2 Subsetting the profileplyr object

The profileplyr object can be subset either by sample, or by rows and columns of the matrix for each sample. This is done using the ‘[ ]’ brackets, with the first position being assay matrix rows, the second position being assay matrix columns, and the third position being the entire matrix for each sample (in addition to the rest of the metadata and range data). For example if you wanted to get the first ten rows and columns of each sample matrix:

proplyrObject[1:10,1:10]

## class: profileplyr 
## dim: 10 10 
## metadata(3): info sampleData rowGroupsInUse
## assays(3): hindbrain1_binned liver1_binned kidney1_binned
## rownames: NULL
## rowData names(3): names score dpGroup
## colnames: NULL
## colData names(0):

The more useful subsetting functionality is likely the ability to subset by samples using the third position of the bracket. For example, if you only wanted the first two samples to export to a deepTools matrix or for further downstream analysis, you could simply just subset with the third position:

proplyrObject[ , , 1:2]

## class: profileplyr 
## dim: 763 100 
## metadata(3): info sampleData rowGroupsInUse
## assays(2): hindbrain1_binned liver1_binned
## rownames: NULL
## rowData names(3): names score dpGroup
## colnames: NULL
## colData names(0):

3.3 Changing sample names

It might be helpful to change the sample names to something that is shorter to make labeling of figures clearer. This can be done separately in most of the visualization packages that a profileplyr object can be used in, but you can change the names within in the sampleData it will be changed for all downstream analyses. After importing either a deepTools matrix or a soGGi object, the sample names are stored as the rownames of the sampleData Dataframe.

rownames(sampleData(proplyrObject))

## [1] "hindbrain1_binned" "liver1_binned"     "kidney1_binned"

The names can simply be changed by reassigning the rownames of this Dataframe.

rownames(sampleData(proplyrObject)) <- c("Hindbrain",
                                         "Liver",
                                         "Kidney")

Importantly, this also changes the names of the matricies stored in the assays() section of the profileplyr object, and the names used for any output method will also be changed (e.g. deepTools, EnrichedHeatmap, ggplot)

3.4 Connecting profileplyr functions using the pipe (%>%) operator

Code involving profileplyr objects can be made clearer using the pipe operator from the magrittr package. This is aided by the fact that the profileplyr object is both the output and input of most functions within the package. This is demonstrated in the Quick Start section, as well as some of the more complex visualizations later in this vignette. The user can go all the way from importing a deepTools matrix or soGGi object, through annotation and grouping, and ending with export to either the generateEnrichedHeatmap() (shown below) or export_deepToolsMat() function

4.1 Export to a deepTools matrix

Any profileplyr object can be exported as a matrix formatted for the deepTools ‘plotHeatmap’ function using the export_deepToolsMat() function. This simply requires the name of the profileplyr object and the path to the desired destination of the gzipped matrix. This function will build the metadata for the deepTools matrix based on the metadata in the profileplyr object. Importantly, the range metadata column of profileplyr object (accessed with rowRanges(object)) that is specified in the ‘rowGroupsInUse’ section of metadata for the whole object (found using metadata(object)$rowGroupsInUse) is automatically used to set the grouping of the exported deepTools matrix. The ‘rowGroupsInUse’ can be changed using the groupBy() function (see more details below).

export_deepToolsMat(proplyrObject, "test_output")

4.2 Convert to an EnrichedHeatmap matrix

To generate a heatmap within R using the EnrichedHeatmap package, a special type of matrix of the ‘normalizedMatrix’ class must be used. The convertToEnrichedHeatmapMat() function takes a profileplyr object as the only required input, and then converts the matrices contained within the assays slot of the profileplyr object to a list of ‘normalizedMatrix’ class objects that can be used in the EnrichedHeatmap() function. See the EnrichedHeatmap vignette for detailed examples for how heatmaps can be concatenated together to visualize all of the data in one figure. Importantly, the metadata stored in the rowRanges slot of the profileplyr object can also be utilized to annotate these heatmaps. NOTE: While this function provides the most flexibility for the user to build a custom set of heatmaps, a more useful function for quick visualization of data stored in a profileplyr object is the generateEnrichedHeatmap() function, which builds these heatmaps with one command using the profileplyr object (see the following section).

library(EnrichedHeatmap)
EH_mat <- convertToEnrichedHeatmapMat(proplyrObject)
EH_mat[[1]]

EnrichedHeatmap(EH_mat[[1]], 
                name = names(EH_mat[1]), 
                column_title = names(EH_mat[1])) +
  EnrichedHeatmap(EH_mat[[2]], 
                  name = names(EH_mat[2]), 
                  column_title = names(EH_mat[2])) +
    EnrichedHeatmap(EH_mat[[3]], 
                  name = names(EH_mat[3]), 
                  column_title = names(EH_mat[3]))

# heatmap not included

4.3 Directly generate a customized EnrichedHeatmap with annotated ranges

profileplyr also contains the generateEnrichedHeatmap() function, which takes a profileplyr object and produces an EnrichedHeatmap that can be easily modified to incorporate metadata from the object. So far the profileplyr object used in this vignette is relatively basic (no groups or additional metadata), but this function generates a multipanel heatmap with a variety of arguments that have been tailored to visualizing the profileplyr object. As we explore more functionality within profileplyr, some of these features will be demonstrated.

By default, if the ranges are divided into multiple groups (e.g. multiple BED files were input) the heatmap will automatically be divided and annotated to reflect these groups. Here we have no groups so we can set the ‘include_group_annotation’ argument to be FALSE. If this were set to TRUE, then colored annotation bars would exist to the left of the heatmaps. The maximum value for the y-axis in the line plots at the top of each heatmap are also automatically set based on the highest mean range signal from all of the samples. This can be set manually with the ‘ylim’ argument, or if ylim = NULL, the maximum will be inferred (i.e. be different) for each individual heatmap.

generateEnrichedHeatmap(proplyrObject, include_group_annotation = FALSE)

4.4 Summarizing signal for ggplot or heatmap visualization

Oftentimes it will be important to use some kind of summary statistic for the entire range when trying to understand and visualize differences between samples. This is often going to be done using the summarize() function. The summarize function requires the name of a profileplyr object, the function used to summarize the ranges (e.g. rowMeans or rowMax), and the type of output. Here the basic options will be demonstrated with the example data which contains no groups, however, note how it is used later in the vignette or in the quick start when additional grouping options are discussed and summarizing the data in this way becomes especially useful.

object_sumMat <- summarize(proplyrObject, 
                           fun = rowMeans, 
                           output = "matrix") 
object_sumMat[1:3, ]

##                               Hindbrain      Liver    Kidney
## chr8_33902396_33903372_genes  0.1878635 0.11450865 0.3375368
## chr17_27678695_27679778_genes 0.6568317 0.08386143 0.1087628
## chr17_23546216_23547490_genes 0.5972133 1.05301136 0.6422771

library(pheatmap)
pheatmap(object_sumMat, 
         scale = "row", 
         show_rownames = FALSE)

proplyrObject_long <- summarize(proplyrObject, 
                                fun = rowMeans, 
                                output = "long") 
proplyrObject_long[1:3, ]

##   dpGroup         combined_ranges    Sample    Signal
## 1   genes  chr8_33902396_33903372 Hindbrain 0.1878635
## 2   genes chr17_27678695_27679778 Hindbrain 0.6568317
## 3   genes chr17_23546216_23547490 Hindbrain 0.5972133

ggplot(proplyrObject_long, aes(x = Sample, y = log(Signal))) + 
       geom_boxplot()

proplyrObject_summ <- summarize(proplyrObject, 
                                fun = rowMeans, 
                                output = "object")

assays(proplyrObject_summ)[[1]][1:3, ]

##      Hindbrain      Liver    Kidney
## [1,] 0.1878635 0.11450865 0.3375368
## [2,] 0.6568317 0.08386143 0.1087628
## [3,] 0.5972133 1.05301136 0.6422771

rowRanges(proplyrObject_summ)[1:3]

## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames            ranges strand |    names     score  dpGroup
##          <Rle>         <IRanges>  <Rle> | <factor> <numeric> <factor>
##   [1]     chr8 33902396-33903372      * |        .         0    genes
##   [2]    chr17 27678695-27679778      * |     ._r1         0    genes
##   [3]    chr17 23546216-23547490      * |     ._r2         0    genes
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths

5.1 K-means and hierarchical clustering of the genomic ranges

Clustering of the signal across the genomic intervals of interest can shed light on patterns that aren’t immediately apparent when looking at the data as a whole. The clusterRanges() function provides a framework to cluster ranges that are contained in a profileplyr object. It should be noted that clustering can also be performed outside of R within the deepTools ‘plotHeatmap’ function, however, profileplyr clustering allows for integration of this clustering with further grouping mechanisms and visualizations within R and profileplyr.

clusterRanges() takes a profileplyr object as its first argument as well as a function to summarize the signal in each range (similar to the summarize() function above). The pheatmap package is then used to cluster the ranges, and the type of clustering depends on whether the user inputs a value for ‘kmeans_k’ (for kmeans clustering) or ‘cutree_rows’ (for hierarchical clustering using hclust). An integer value entered for either of these arguments will specify the number of clusters used for each method. If both ‘kmeans_k’ or ‘cutree_rows’ are left as NULL (default), then a heatmap will be printed with hierarchical clustering, but no distinct clusters defined, and no profileplyr object will be returned. This might be helpful as an initial and quick look at the ranges or as a means to determine the number of clusters to try.

clusterRanges(proplyrObject, 
              fun = rowMeans)

# this code prints heatmap (does not return), but heatmap not shown here to save space

If an integer is entered for either ‘kmeans_k’ or ‘cutree_rows’, then a profileplyr object will be returned with a cluster assigned to each range, and a column is added to the range metadata with the cluster for each range. Further, the ‘rowGroupsInUse’ is changed to this column, meaning that if this object is exported as a deepTools matrix with the export_deepToolsMat() function, the heatmap generated by ‘plotHeatmap’ will be grouped by the clusters.

set.seed(0)
kmeans <- clusterRanges(proplyrObject, 
                        fun = rowMeans, 
                        kmeans_k = 4)

## [1] "K means clustering used."
## [1] "A column has been added to the range metadata with the column name 'cluster', and the 'rowGroupsInUse' has been set to this column."

rowRanges(kmeans)[1:3]

## GRanges object with 3 ranges and 4 metadata columns:
##       seqnames            ranges strand |    names     score  dpGroup
##          <Rle>         <IRanges>  <Rle> | <factor> <numeric> <factor>
##   [1]     chr8 33902396-33903372      * |        .         0    genes
##   [2]    chr17 27678695-27679778      * |     ._r1         0    genes
##   [3]    chr17 23546216-23547490      * |     ._r2         0    genes
##         cluster
##       <ordered>
##   [1]         4
##   [2]         3
##   [3]         1
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths

metadata(kmeans)$rowGroupsInUse

## [1] "cluster"

To visually inspect the clusters, the heatmap can also be printed (though it will not be returned) if the ‘silent’ argument is set to FALSE. A profileplyr object will still be returned even if silent is set to FALSE.

Kmeans clustering:

set.seed(0)
kmeans <- clusterRanges(proplyrObject, 
                        fun = rowMeans, 
                        kmeans_k = 4, 
                        silent = FALSE)

## [1] "K means clustering used."

## [1] "A column has been added to the range metadata with the column name 'cluster', and the 'rowGroupsInUse' has been set to this column."

Hierarchical Clustering:

hclust <- clusterRanges(proplyrObject, 
                        fun = rowMeans, 
                        cutree_rows = 4, 
                        silent = FALSE)

## [1] "Hierarchical clustering used. It is advised to avoid this option with large matrices as the clustering can take a long time. Kmeans is more suitable for large matrices."

## [1] "A column has been added to the range metadata with the column name 'cluster', and the 'rowGroupsInUse' has been set to this column."

5.2 Output matrix with cluster information to deepTools

The profileplyr object that is the output of clusterRanges() can be exported as a deepTools matrix, which can then be used as an input for plotHeatmap.

export_deepToolsMat(kmeans, "kmeans_cluster_mat")

This matrix can then be passed directly to ‘plotHeatmap’, either by using the terminal or by using the system() function in R.

Code from command line :

plotHeatmap -m kmeans_cluster_mat.gz -o kmeans_cluster_mat.jpg –samplesLabel Hindbrain Liver Kidney –startLabel start –endLabel end –xAxisLabel “distance (bp)”

5.3 Generate group-annotated heatmap in R directly with generateEnrichedHeatmap()

After clustering the profileplyr object can be passed directly as an argument into the generateEnrichedHeatmap() function, and by default the heatmap will be grouped and annotated by these clusters, which were automatically set as the ‘rowGroupsInUse’ in the clusterRanges() function. Assuming that the ‘include_group_annotation’ argument of this function is set to the default value of TRUE, whichever metadata column is set to the ‘rowGroupsInUse’ slot will be used for the grouping and annotation seen below.

generateEnrichedHeatmap(kmeans)

5.4 Visualize mean range signal for each cluster with ggplot

An alternative way to visualize the clusters is by summarizing the signal over each range and then plotting those means by cluster. The flexibility with grouping and annotating within profileplyr makes this relatively easy. This is also a good demonstration of using piped code, going all the way from import to clustering, then to summarizing by mean range signal, and finishing with ggplot visualization:

library(magrittr)

set.seed(0)
import_deepToolsMat(example) %>%
  clusterRanges(fun = rowMeans, 
                kmeans_k = 4, 
                silent = TRUE) %>% 
  summarize(fun = rowMeans, 
            output = "long") %>% 
  ggplot(aes(x = Sample, y = log(Signal))) + 
         geom_boxplot() + 
         facet_grid(~cluster) + 
         theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
         scale_x_discrete(labels= c("Hindbrain", "Liver", "Kidney"))

5.5 Gene annotation of ranges

Understanding the genes that are in the proximity to a particular set of ranges allows for potential understanding of functional consequences of signal patterns within these ranges. By annotating the ranges within a profileplyr object with genes, the user can connect signal of particular samples to specific genes and the functions attributed to those genes. profileplyr provides two methods of annotating the ranges within the object using two established Bioconductor packages, ChIPseeker and rGREAT.

anno <- annotateRanges(proplyrObject, 
                       TxDb = "mm10")

## >> preparing features information...      2019-03-29 01:03:18 
## >> identifying nearest features...        2019-03-29 01:03:19 
## >> calculating distance from peak to TSS...   2019-03-29 01:03:19 
## >> assigning genomic annotation...        2019-03-29 01:03:19 
## >> adding gene annotation...          2019-03-29 01:03:28 
## >> assigning chromosome lengths           2019-03-29 01:03:29 
## >> done...                    2019-03-29 01:03:29

rowRanges(anno)[1:3]

## GRanges object with 3 ranges and 16 metadata columns:
##       seqnames            ranges strand |    names     score  dpGroup
##          <Rle>         <IRanges>  <Rle> | <factor> <numeric> <factor>
##   [1]     chr8 33902396-33903372      * |        .         0    genes
##   [2]    chr17 27678695-27679778      * |     ._r1         0    genes
##   [3]    chr17 23546216-23547490      * |     ._r2         0    genes
##                                      annotation   geneChr geneStart
##                                     <character> <integer> <integer>
##   [1]  Intron (uc009lkk.2/19663, intron 1 of 8)         8  33782644
##   [2]  Intron (uc008bpj.2/23969, intron 1 of 9)        17  27685239
##   [3] Exon (uc033hbl.1/uc033hbl.1, exon 1 of 1)        17  23550799
##         geneEnd geneLength geneStrand      geneId transcriptId distanceToTSS
##       <integer>  <integer>  <integer> <character>  <character>     <numeric>
##   [1]  33929863     147220          2       19663   uc009lkk.2         26491
##   [2]  27711117      25879          1       23969   uc008bpm.2         -5461
##   [3]  23561634      10836          1      320020   uc056zdm.1         -3309
##                  ENSEMBL        SYMBOL
##              <character>   <character>
##   [1] ENSMUSG00000031586         Rbpms
##   [2] ENSMUSG00000040276       Pacsin1
##   [3] ENSMUSG00000117322 6330415G19Rik
##                                                        GENENAME
##                                                     <character>
##   [1]           RNA binding protein gene with multiple splicing
##   [2] protein kinase C and casein kinase substrate in neurons 1
##   [3]                                RIKEN cDNA 6330415G19 gene
##       annotation_short
##              <ordered>
##   [1]           Intron
##   [2]           Intron
##   [3]             Exon
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths

generateEnrichedHeatmap(anno, 
                        extra_annotation_columns = "annotation_short",
                        include_group_annotation = FALSE)

anno_great <- annotateRanges_great(proplyrObject, 
                                   species = "mm10")
rowRanges(anno_great)[1:3]

## GRanges object with 3 ranges and 5 metadata columns:
##       seqnames            ranges strand |    names     score  dpGroup
##          <Rle>         <IRanges>  <Rle> | <factor> <numeric> <factor>
##   [1]     chr1 14181321-14182231      * |   ._r487         0    genes
##   [2]     chr1 14181321-14182231      * |   ._r487         0    genes
##   [3]     chr1 16201828-16202701      * |   ._r287         0    genes
##            SYMBOL distanceToTSS
##       <character>     <numeric>
##   [1]        Eya1        128423
##   [2]        Xkr9        513005
##   [3]       Rdh10         96491
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths

6.1 Switch output grouping columns within existing range metadata

If the ‘group’ argument is a character string then it must match a name of one of the columns in mcols(proplyrObject). In this case the groupBy function will change the column that will be used for grouping a profileplyr object during visualization. While many of the profileplyr functions shown above will set the grouping column (specified in the ‘rowGroupsInUse’ section of metadata(proplyrObject)) to the appropriate column for that function, the user may want to use another column for the grouping of the deepTools matrix to be output. The user may have added an additional column to the range metadata that can be used, or a column that is generated in one of the above functions, but is not the default grouping column, might be useful. For example, the default grouping for the annotateRanges() function is a combination of the inherited groups and the annotations, however, maybe the user wants to disregard the inherited grouping of ranges and group all of the ranges by the new annotation. In this case they would change the grouping to the ‘annotation_short’ column, as shown below.

anno <- annotateRanges(proplyrObject,
                       TxDb = "mm10")

metadata(anno)$rowGroupsInUse

## [1] "dpGroup"

anno_switchGroup <- groupBy(anno, 
                            group = "annotation_short")
metadata(anno_switchGroup)$rowGroupsInUse

## [1] "annotation_short"

6.2 Group by user-supplied GRanges

If the ‘group’ argument is a GRanges or a GRangesList then the overlap between each GRanges in this list and the ranges of the profileplyr object will be determined. A column will be added to the range metadata (column name is ‘group_and_overlap’) that shows a combination of the inherited group and whether the range overlaps with any of the GRanges. The ‘rowGroupsInUse’ will be changed to this column by default. If the inherit_groups argument is set to FALSE, then only the overlaps with the GRanges will determine the grouping.

One consideration with this function is what should be done with ranges that overlap multiple GRanges. The default behavior of this function is to duplicate the range and place it in each overlap group for each GRanges. If separateDuplicated = TRUE, then separate groups will be generated for the ranges that overlap multiple GRanges so that every range is in the resulting profileplyr object exactly once. This option is discouraged with high numbers of input GRanges, as the combinations of overlapping groups can get quite large. Lastly, the non-overlapping regions can be included in the heatmap if the ‘include_nonoverlapping’ argument is set to TRUE.

Here we import BED files containing the top 5000 H3K27ac peaks from hindbrain and liver using the rtracklayer package. These peaks sets were adapted from BED files that were downloaded from ENCODE using the provided links. The GRanges are then combined into a GRangesList to be used in the ‘group’ argument for the groupBy() function.

library(rtracklayer)

# import bed files as GRanges
hindbrain_K27ac <- import.bed(system.file("extdata", 
                                          "K27ac_hindbrain_top5000.bed", 
                                          package = "profileplyr") )
liver_K27ac <- import.bed(system.file("extdata", 
                                      "K27ac_liver_top5000.bed", 
                                      package = "profileplyr") )
# make GRangesList
K27ac_GRlist <- GRangesList(hindbrain_K27ac, 
                            liver_K27ac)
K27ac_GRlist_names <- c("hindbrain_K27ac", 
                        "liver_K27ac")

# subset the profileplyr object to just hindbrain and liver, 
# and group by GRanges
K27ac_groupByGR <- proplyrObject[,,grepl("Hindbrain|Liver", 
                                         names(assays(proplyrObject)))] %>%
                   groupBy(group = K27ac_GRlist,
                           GRanges_names = K27ac_GRlist_names, 
                           include_nonoverlapping = FALSE,
                           inherit_groups = FALSE)

rowRanges(K27ac_groupByGR)[1:3]

## GRanges object with 3 ranges and 8 metadata columns:
##       seqnames            ranges strand |    names     score  dpGroup
##          <Rle>         <IRanges>  <Rle> | <factor> <numeric> <factor>
##   [1]    chr16 23520075-23521344      * |     ._r4         0    genes
##   [2]    chr17 55445046-55446488      * |    ._r15         0    genes
##   [3]    chr16 77399830-77401368      * |    ._r21         0    genes
##       overlap_matrix GRoverlap_names   group_and_GRoverlap        name
##             <matrix>       <ordered>             <ordered> <character>
##   [1]            1:0 hindbrain_K27ac genes_hindbrain_K27ac        <NA>
##   [2]            1:0 hindbrain_K27ac genes_hindbrain_K27ac        <NA>
##   [3]            1:0 hindbrain_K27ac genes_hindbrain_K27ac        <NA>
##         score.1
##       <numeric>
##   [1]       132
##   [2]       145
##   [3]       194
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths

As with any profileplyr object, this can be exported to deepTools ‘plotHeatmap’ using the profileplyr function export_deepToolsMat()

export_deepToolsMat(K27ac_groupByGR, 
                    con = "K27ac_GRoverlap.gz")

This matrix can then be passed directly to ‘plotHeatmap’, either by using the terminal or by using the system() function in R.

Code from command line:

plotHeatmap -m K27ac_GRoverlap.gz -o K27ac_GRoverlap.jpg –samplesLabel Hindbrain Liver Kidney –startLabel start –endLabel end –xAxisLabel “distance (bp)”

The profileplyr object can also be passed directly to the generateEnrichedHeatmap() function to be visualized in R.

generateEnrichedHeatmap(K27ac_groupByGR)

6.3 Group by user-supplied gene list

If the ‘group’ argument is a list (i.e. class(group) returns “list”), then it is assumed that this list contains lists of gene sets. These gene sets within this list can either be character vectors of gene symbols, or data frames with the gene symbols as rownames. Importantly, if all of the elements of the list are data frames AND they all have the same columns (as determined by having matching column names), then those columns will be added to the range metadata to potentially be used to annotate the ranges. This is especially useful for annotating the ranges with additional gene-level metadata (e.g. gene expression). However, if any of the elements of the list are character vectors, or if they are data frames with columns that are different, then the additional columns will not be included in the range metadata and only the overlap information will be stored.

The gene symbols will be used to determine overlap with the genes that are associated with each range of the profileplyr object. This type of grouping will typically follow one of the annotation functions: annotateRanges() or annotateRanges_great(). The names of the list elements (names(list)) will be used as the names of each set of genes in the range metadata and any exported deepTools matrix. A new column will be added to the range metadata that shows a combination of the inherited group and whether the range overlaps with any of the gene sets.

gene_list_hindbrain_liv <- system.file("extdata", 
                                            "gene_list_hindbrain_liv.RData", 
                                            package = "profileplyr") 

gene_list_character <- readRDS(gene_list_hindbrain_liv) 
names(gene_list_character)

## [1] "upInHindbrain_VS_Liver"   "downInHindbrain_VS_Liver"

head(gene_list_character[[1]])

## [1] "Tuba1a" "Agap1"  "Nnat"   "Col2a1" "Rtn1"   "Map1b"

proplyrObject_HB_Liv <- proplyrObject[ , , grepl("Hindbrain|Liver", 
                                                 names(assays(proplyrObject)))]

proplyrObject_HB_Liv <- annotateRanges_great(proplyrObject_HB_Liv, 
                                             species = "mm10") 

proplyrObject_HB_Liv <- groupBy(proplyrObject_HB_Liv, 
                                group = gene_list_character)

rowRanges(proplyrObject_HB_Liv)[1:3]

## GRanges object with 3 ranges and 8 metadata columns:
##       seqnames            ranges strand |    names     score  dpGroup
##          <Rle>         <IRanges>  <Rle> | <factor> <numeric> <factor>
##   [1]     chr1 16201828-16202701      * |   ._r287         0    genes
##   [2]     chr1 20867186-20868863      * |   ._r211         0    genes
##   [3]     chr1 24107322-24109025      * |    ._r88         0    genes
##            SYMBOL distanceToTSS overlap_matrix        GLoverlap_names
##       <character>     <numeric>       <matrix>              <ordered>
##   [1]       Stau2        317041            1:0 upInHindbrain_VS_Liver
##   [2]       Paqr8        -22597            1:0 upInHindbrain_VS_Liver
##   [3]      Col9a1        -69436            1:0 upInHindbrain_VS_Liver
##                group_and_GLoverlap
##                          <ordered>
##   [1] genes_upInHindbrain_VS_Liver
##   [2] genes_upInHindbrain_VS_Liver
##   [3] genes_upInHindbrain_VS_Liver
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths

export_deepToolsMat(proplyrObject_HB_Liv, 
                    con = "great_groupByGene_HB_Liv")

gene_list_hindVliver_df <- system.file("extdata", 
                                       "gene_list_hindbrain_liv_df.RData", 
                                       package = "profileplyr") 

gene_list_dataframe <- readRDS(gene_list_hindVliver_df) 
names(gene_list_dataframe)

## [1] "upInHindbrain_VS_Liver"   "downInHindbrain_VS_Liver"

head(gene_list_dataframe[[1]])

##        hindbrain_liv_stat
## Tuba1a           75.80852
## Agap1            53.32196
## Nnat             51.68476
## Col2a1           49.94588
## Rtn1             48.51060
## Map1b            44.71733

example <- system.file("extdata", 
                       "example_deepTools_MAT", 
                       package = "profileplyr")

# NOTE: could also use BamBigwig_to_chipProfile() and as_profileplyr()
# within R to start from bigwig/BAM and BED file. 

import_deepToolsMat(example) %>%
  .[ , , grepl("hindbrain|liver", names(assays(.)))] %>%
  annotateRanges_great(species = "mm10") %>%
  groupBy(group = gene_list_dataframe,
          inherit_groups = FALSE) %>%
  generateEnrichedHeatmap(extra_annotation_columns="hindbrain_liv_stat",
                          sample_names = c("Hindbrain", 
                                           "Liver"),
                          extra_anno_width = 8)